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ABSTRACT 

We provide a both qualitative and quantitative comparison among different 
approaches aimed to solve the problem of non-linear diffusive acceleration of 
particles at shocks. In particular, we show that state-of-the-art models (nu- 
merical, Monte Carlo and semi-analytical) , even if based on different physical 
assumptions and implementations, for typical environmental parameters lead 
to very consistent results in terms of shock hydrodynamics, cosmic ray spec- 
trum and also escaping flux spectrum and anisotropy. Strong points and limits 
of each approach are also discussed, as a function of the problem one wants 
to study. 
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1 FROM DSA TO NLDSA 



In th e late 70 's m any authors (JKriniskiil 



1977 



Axford. Leer fc Skardon 



1977 



Blandford fc Ostriker 



1978 



Bell 



19781) introduced the test-particle theory of particle acceleration at strong coUi- 



sionless shocks due to first order Fermi mechanism. However, quantitative estimates soon 
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2 Caprioli, Kang, Vladimirov and Jones 

pointed out that this mechanisin, called Diffusive Shock Acceleration (DSA), may be so 
efficient that the back-reaction of the accelerated particles on the shock dynamics can not 
be neglected. 

The obvious theoretical challenge is how to model effectively the full shock dynam- 
ics. Modeling shocks with efficient particle acceleration using plasma simulations (PIC and 
hybrid) is extremely computationally expensive for several reasons. First, the energies of 
charged particles participating in the process range from the low thermal energies of cold 
plasma to the ultra-relativistic energies of cosmic rays, and bot h extremes of particle spectra 
are d ynamically important if acceleration is efficient (see, e.g.. I Vladimirov. Ellison &: Bykov 
20081 . for estimates of computational requirements). Second, simulations need to be done in 
three dimensions because of the possibility of non-physical suppression of importarit pro - 
cesses in one- and two-dimensional simulations (see, e.g., iJones. Jokipii &: Baring! Il998[ l. 
Therefore, approximate methods must be used to model efficient particle accelerating shocks. 

The first step towards a Non-Linear theory of DSA (NLDSA), in which particle accel- 
eration and shock dynamics are calcula t ed se lf-consistently, is represented by the two-fluid 



model proposed by 



Drury fc Volkl (jlQSOl . Il98ll ) . This approach is indeed very useful to study 



the dynamics of a shock modified by the presence of accelerated particles (treated as a diffu- 
sive, relativistic fluid) but by design it cannot provide any information about the spectrum 
of the cosmic rays (CRs) which are produced. 

A kinetic study of the NLDSA at astrophysical shocks also requires a statistical descrip- 
tion of transport of superthermal particles. This relevant piece of Physics has been taken 
into account by adopting a number of rather different approaches; namely 



• fully numerical simulations, in which a time-dependent diffusion- convection equation for 
the CR transport is solved in conce rt with coupled gas dy n amics conservation laws , like i n 



Belli (Il987h: 



Berezhko fc Vo 



Kang fc JonesI (12005 



si J2005I . bood . 



klfll997h: 



Kangfc JonesI (Il997h: 



20071 ): IZirakashvih fc Aharonianl fl201 



Kang. Jones fc Gieseleii (l2002h : 



m 



stationary Monte Carlo numerical simulations of the full par ti cle population, like in 



Ellison fc Eichlei ( 



1984 



1985) 



Ellison. Mobius &: Paschmann (11990) 



Ellison. Baring fc Jones 



fll995l . Il996h 



Ellison fc Double! ( 120021 ): 



Vladimirov. Ellison fc Bykov! ( l2006l ) 



semi-analytic solutions of the stationary (or quasi- stationary, see lBlasi. Amato fc Caprioli 



20071 ) diffusion-convection equation coupled to the gas dynamic equations, like in 



Malkov 
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1997); 



Malkov. Diamond fc Volkl fl2000h: 



Caprioli et al. 



(J2nn9ah 



Blasi 



Caprioli. Amato. Blasi 



(2002 



feoiObI ). 



20041 ): lAmato fc Blasil fcoosl . l2006f l: 



All these approaches make consistent predictions about the main consequences of the 
shock modification, i.e. the appearance of an upstream precursor created by the CR pressure 
around the shock, which slows down the incoming fluid. As a consequence of this CR-induced 
velocity gradient, the fluid compression ratio is no longer 4 for any strong shock. Hence, 
the test-particle prediction that the spectrum of accelerated particles has to be oc p~^ is 
no longer valid. Nevertheless, it is still true that the spectral slope mainly depends on the 
compression factor actually "felt" by CRs. Since the larger is its momentum, the farther from 
the shock a particle can diffuse, high (low) momentum particles probe a total compression 
ratio larger (smaller) than 4. The resulting CR spectrum becomes rather concave, being 



harder (softer) than p '^ at the higher 



shocks the reader can refer to 



(ll99ll); 



Malkov fc Drurvl f l200ll ) 



■ower) energies. For accurate review s of CR-modifled 



Drurvl fll983f ): iBlandford fc Eichleil fll987h : 



Jones &: Ellison 



In this paper we want, for the flrst time, to construct both qualitative and quantita- 
tive comparisons among different approaches to NLDSA, in order to highlight their analo- 
gies and differences. In particular, we compare three state-of-the-art kinetic methods which 
are representa t ive of the three categories listed abov e: the fully numerical calculat i on by 
Kang fc Jones! (120071 ). the Monte C arlo simulation bv I Vladimirov. Ellison fc Bykovl (120061 ) 



and the semi-analytical solution by 



Caprioli. Amato. Blasi 



(l2010b[ ). The flnal goal of this 



work is to provide a useful tool illustrating which physical aspects of the problem can be 
taken into account within each approach, which assumptions are a posteriori justifled and 
consistent with any observational signature, and flnally what is the computational effort 
required to tackle a given problem. 

The present work is organized as follows: the next three sections are dedicated to the 
description of the three different techniques adopted to solve the problem of NLDSA. For 
each of them, the assumptions (e.g. time-dependence, isotropy of the CR distribution func- 
tion) and the equations of the model, the injection of particles in the acceleration mechanism 
and the turbulent heating in the precursor are outlined. In ^ the results obtained within 
different approaches are compared, in terms of both hydrodynamics and CR distribution 
function. We comment and conclude in ^ 
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4 Caprioli, Kang, Vladimirov and Jones 
2 NUMERICAL APPROACH 

2.1 Basic equations 

In their kinetic simulations of DSA, Kang and Jones (hereafter KJ) evolve the standard 
time-dependent gasdynamic conservation laws with the CR pressure terms in Eulerian form 



for o ne- dimensional plane-parallel geometry (JKang. Jones fc Gieseler 
2003), 



2002 



Kang fc Jones 



d_ 
di 



dp 

m 

(pu) 



-^{P^' + P^ + Pc) 



d_ 

di 



1 2 Pg 

2 7- 



d 



JUPg 

7-1 



u- 



dPc 

dx 



+ >V(x,t) - C{x,t) . 



(1) 
(2) 
(3) 



where p = p{x, t) is the gas density, u = u{x, t) is the fluid velocity, P = P{x, t) is the gas (Pg) 
or CR {Pc) pressure (eq. |22]) and 7 = 5/3 is the gas adiabatic index. C{x,t) is the injection 
energy loss term, which accounts for the energy carried away by the suprathermal particles 
injected into the CR component at the subshock and is subtracted from the postshock gas 
immediately behind the subshock. The gas heating due to the A lfven wave dissipation in 
the upstream region is modeled following iMcKenzie fc Volkl (119821 ) and is given by the term 

dPr 



W{x,t) = -VA{x,t)- 



dx 



(4) 



where VA{x,t) = Bo/^yAnp{x,t) is the local Alfven speed and Bq is the magnetic field 
strength far upstream. Here and in the following we will label with the subscript quantities 
at upstream infinity, and with 1 and 2 quantities immediately upstream and downstream of 
the subshock, respectively. These equations can be used to describe parallel shocks, where 
the large-scale magnetic field is aligned with the shock normal and the pressure contribution 
from the turbulent magnetic fields can be neglected. 

The CR population is evolved by solving the diffusion-convection equation for the pitch- 
angle-averaged (isotropic in momentum space) distribution function, f{x,p,t), in the form: 



dl 
dt 



+ [u + u^ 



dx 



d_ 
dx 



D{x,p) 



dl 
dx 



+ 



p df d{u + Uw) 
3 dp dx 



Skilling 



1975 



(5) 



for a derivation) and 



ly rewritten and solved 



where D(x,p) is the spatial diffusion coefficient (see e.g 

p is the scalar, momentum magnitude. This equation is then usefu 

i n ter ms of g{x,p,t) = p'^f^XjPjt) and y = log(p) as explained in 

( |2002| ). Eq. 7. The CR population is isotropized with respect to the local Alfvenic wave 

turbulence, which would in general move at a speed u^j with respect to the plasma. Since 



Kang. Jones fc Gieseler 
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Comparison of Different Methods for NLDSA 5 

the Alfven waves upstream of the subshock are expected to be estabhshed by the streaming 
instabihty, the wave speed is set there to be u^ = —va- Downstream, it is hkely that the 
Alfvenic turbulence is nearly isotropic, hence Uy^ = there. 

K J typically initialize their simulations by creating pure hydrodynamical shocks and then 
evolving them along with CRs. The shocks described here were initialized from Rankine- 
Hugoniot (R-H) relations, although alternatives, such as reflecting boundaries can be em- 
ployed. It is important to note that the initial shock speed established by the R-H conditions 
(and thus, also the Mach number) will change as the shock evolves. That results from in- 
creased compression through the shock. This feature, which is not a factor in the two steady 
state shock models, complicates somewhat our comparison methodology. 

Previous studies showed that CR modified shocks evolve to a dynamical self-similar 
state, with constant total compression and downstream pressure balance; however, the 
CR population continues to stretch to ever higher momenta unless high energy CRs es- 
cape either through an up per momentum boundary, or a spatial "free escape" boundary 
( iKang. Ryu fc Jonesll2009l ). The former is simply implemented by setting f{x,p) = for 
P > Pmax, while the latter (employed in this study) is effected by setting f{x,p) = for 
X > Xq. Then the modified shocks and their CR population will evolve until CR escape 
through either of these boundaries balances DSA within the shock structure. From that 
time forward the shock and its CR population become steady and the spectrum of escaping 
particles can be worked out numerically as 
'df 



(f)escip) = -D{p) 



dx 



(6) 

XO 



2.2 Method of solution 

Equations [H [2], [31 and E] are simultaneously integrated by finite difference methods in the 
CRASH (Cosmic-Ra y Acceleration SHock) c o de. T he detailed description of the CRASH 
code can be found in iKang. Jones fc Gieselerl ( 20021 ). Three features of CRASH are impor- 



tant to our discussion below. First, CRASH applies an adaptive mesh refinement technique 
around the subshock. So the precursor structure is adequately resolved to couple the gas to 
the CRs of low momenta, whose diffusion lengths can be at least several orders of magnitude 
smaller than the precursor width. Second, CRASH uses a subgrid shock tracking; that is, the 
subshock is a true discontinuity, whose position is followed accurately within a single cell on 
the finest mesh refinement level. Consequently, the effective numerical subshock thickness 
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6 Caprioli, Kang, Vladimirov and Jones 

needed to compute the spatial derivatives in Eq. [5] is always less than the single cell size of 
the finest grid. Third, the exact subshock speed is calculated at each time step to adjust the 
rest frame of the simulation, so that the subshock is kept inside the same grid cell through- 
out. These three features allow one to obtain good numerical convergence in the solutions 
with a minimum of computational effort. A typical shock simulation with CRASH evolved 
to steady state with CRs accelerated to ~ lO'^GeV/c requires roughly 300 CPU hours on 
current microprocessors. CRASH does not currently have a parallel implementation. The 
evolving shock simulation time can be shortened by as much as an or der of magnitud e using 
the finite volume, CGMV algorithm for CR transport discussed in iJones fc KangI ( l2005l ) 
instead of traditional finite difference methods. However, thermal leakage in the sophisti- 
cated form discussed in the following subsection has not been implemented with the CGMV 
algorithm. So, since one of our objectives here is comparisons of injection methods, we do 
not employ CGMV. 



2.3 Thermal leakage injection 

Since the velocity distribution of suprathermal particles in the Maxwellian tail is not isotropic 
in the shock frame, the diffusion-convection equation cannot directly follow the injection 
from the non-diffusive thermal pool into the diffusive CR population. In the CRASH code, 
thermal leakage injection is emulated numerically by adopting a transparency function, 
Ti,sc{eB,v), which expresses the probability for downstream thermal particles at given ran- 
dom velocity, v, to swim successfully upstream across the subshock through the postshock 
cyclotron- instability-generate d MHD waves, whose amplitude is parameterized by es (see 



Kang. Jones &: Gieselerl 120021 . Eq. 8). The condition of finite probability for suprathermal 
downstream particles to cross the subshock (i.e. Tcsc > for p larger than an injection 
momentum pinj ~ (1 + 1.07/eB)'mpU2 effectively selects the lowest momentum of the parti- 
cles entering the CR population. The model transparency function increases smoothly from 
zero to unity over pi^j < p <~ 5pinj. Consequently, the transition from the Maxwellian to 
the quasi-power-law distribution develops smoothly in momentum space. The parameter, 
es = Bo/B±, is the ratio of the magnitude of the large-scale magnetic field aligned with 
the shock normal, Bq, to the amplitude of the postshock wave field that interacts with low 



energy particles, B±. This is an adjustable parameter, although, as discussed in 



Malkov 



(1l997[ ). it is confined theoretically and experimentally to values, e^ ~ 0.3. 
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Comparison of Different Methods for NLDSA 7 
3 MONTE CARLO APPROACH 

3.1 Model assumptions 

The method of solving the NLDSA problem based on the Monte Carlo particle transport, 
developed by Ellison and co-workers and later adapted for modeling magnetic field amplifi- 
cation in shocks by Vladimirov (hereafter EV), makes the following physical assumptions: 

• all quantities depend only on one spatial coordinate, but the flow velocity and/or the 
uniform magnetic field may have a component along the perpendicular direction; 

• a steady state is realized; 

• mass, momentum and energy conservation laws determine the nonlinear shock structure; 

• the transport of all particles may be described statistically. In the simplest case, it 
is assumed that the momentum vector of every particle gets randomly perturbed in the 
course of propagation in such a way that particle distribution is isotropized in the course 
of a momentum-dependent collision time, which is a parameter of the model. If scattering 
rate is constant in space and time, then on large enough scales such particle transport 
models diffusion. Gyration around a uniform magnetic field may be included along with the 
stochastic component of the motion; 

• particles of all energies are treated on an equal basis, i.e. there is no need for a strict 
distinction between thermal and superthermal particles. However, such a distinction may be 
artificially imposed (see §3.3p : 

• collisions that facilitate particle diffusion are elastic in the reference frame moving at a 
speed Uw with respect the thermal plasma, where u^ is chosen the same way as in the other 
models. 

The problem of NLDSA is thus stated as that of finding a steady state particle distri- 
bution function in a system where a supersonic motion exists at an infinite boundary, and 
the fiow is subsonic at the opposite boundary. The distribution function F{x, p) in this 
case retains the information about the direction of the momentum vector p and spans the 
range of momenta from thermal to Pmax- Here and in the following we indicate with F{x, p) 
the whole distribution function (thermal -|- cosmic rays) and with /(x,p) the cosmic ray 
distribution only. 

© 2010 RAS, MNRAS OOO.ITIMI 



8 Caprioli, Kang, Vladimirov and Jones 
3.2 Method of solution 

The method by which this problem is solved involves an iterative application of the following 
two steps: 

1) a trial flow speed u{x) is chosen based on the conservation laws and on the best-yet 
estimate of F(x, p), and 

2) the updated particle distribution F(a;, p) is calculated with the Monte Carlo technique 
based on the trial flow speed u{x). 

These steps are repeated until mass, momentum and energy conservation laws hold for the 
particle distribution F(a;, p) calculated with the chosen u{x). 

In practice, the flrst step is realized by adjusting the previous guess for u{x) using the 
conservation laws. The conservation of momentum is used to choose u{x) in the supersonic 
region (i.e., upstream), and the conservation of momentum and energy is used to choose 
the flow speed in the subsonic region (i.e. downstream), taken as uniform. The equations 
expressing the conservation of mass, momentum and energy fluxes can be written in terms 
of -F(x, p) as 



mpV^F{x, p) d p = pqUq (7) 

PxVxF{x, p)d^p = poul + Pgfi , (8) 

/I . 7 
K{p)v^F{x, p)d^p = -pou'o H rUoPgfi - Qcsc , (9) 

2 7-1 

where the subscript x labels the component along the axis x, K{p) is the kinetic energy of 

a particle with momentum p and Qesc is the energy flux carried away by particles escaping 

the system (see below for details). This way of dealing with a global distribution function is 

in contrast to writing the left hand side as the sum of gas and CR pressures, which would 

require F{x, p) to be isotropic in p-space. 

The second step is realized by introducing a population of thermal particles either far 

upstream, or at a location xa (deflned below), and propagating them in time steps small 

compared to the collision time, until they escape. A particle is considered to have escaped 

if a) it is found many diffusion lengths downstream, or b) it crosses the upstream free 

escape boundary. As the particles are propagated, their momentum are randomly perturbed 

as described in the next paragraph. When particles cross selected coordinates Xj, their 

momenta and statistical weights are summed to calculate the particle distribution function 

in a momentum bin: 

© 2010 RAS, MNRAS OOO.nH24l 
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peApki^eAfMi -f^*: ^'^ ^'■ 
Here, the suinination takes place over the events of a particle crossing the coordinate Xj: 
Apk is a bin in momentum space, Afii is a bin in angles (/i = Px/p) and cij is the properly 
normalized statistical weight of the particle. In a similar fashion, the code calculates the 
moments of the distribution function, including the momentum and energy fluxes defined in 

Eq.[8]and[i 

When particles are propagated, they move in straight lines for a time At equal to a frac- 
tion of the collision time t^. After that, the momentum in the local plasma frame is rotated 
by a random angle chosen so that the averaged particle motion represents diffusion with the 
appropriate diffusion coefficient. The local plasma frame for CR particles is assumed to be 
moving at a speed u{x) — va{x) with respect to the subshock, where va{x) = Bq/ ^jATip{x) 
is the local Alfven speed. For thermal particles, the frame in which momentum rotation 
is performed is moving at a speed u{x ) with respect to the subshock. The details of this 
procedure are described in IVladimirovl (120091 ). In regions where the particle distribution is 
isotropic (e.g., far downstream), the code may set At = tc, which simulates a random walk 
process and corresponds to the diffusion approximation. 

A gradient of the local plasma speed in the above procedure leads to the first order Fermi 
acceleration of particles. This applies to the thermal particles as well as to the superthermal 
ones. Injection of thermal particles into the acceleration process occurs in this model only 
due to their Fermi-I acceleration at the subshock. 

Elastic scatterings of thermal particles in the shock precursor will, in the presence of a 
small fiow speed gradient, lead to adiabatic heating: Pg oc p"' (this heating is, in a sense, 
due to the same Fermi-I acceleration process, but on a smooth variation of the fiow speed 
instead of a discontinuity). If additional heating operates in the precursor due to some wave 
dissipation processes, the treatment of thermal particles in the precursor is modified as 
described in §3.31 

In the Monte Carlo particle transport routine, the escape of energetic particles upstream 
is modeled by removing from the simulation the particles that cross the coordinate of the 
free escape boundary, xq, in the upstream direction. The value of xq is a free parameter of 
the model. The energy fiux carried away by the escaping particles is refiected in the value of 
the energy fiux downstream of Xq. The amount of escaping energy fiux, Qesc^ is calculated 
in the EV method as the difference between the energy fiux upstream of Xq and the energy 

© 2010 RAS, MNRAS 000.ITH241 



10 Caprioli, Kang, Vladimirov and Jones 

flux downstream of Xq (see equation ([9])), wfiile tlie spectral density of the flux of escaping 
particles can be calculated in the simulation as 

4>esc{p) = dfl CfiF{Xo,p,fi). (11) 

As a final remark, the computational time required by this Monte Carlo approach is on 
the order of tens to hundreds of processor hours, and parallel processing is implemented in 
the latest version. 

3.3 Precursor heating 

In order to include precursor heating, particles are introduced at an upstream location 
xa, which is close to the subshock, but not too close, in order that the thermal particle 
distribution is still isotropic at xa- Upstream of xa, instead of calculating F{x, p) and its 
moments like Pg{x) from particle propagation, the model finds Pg{x) by solving the equation 



p^{x) d 
u[x) 



'J — Idx 



p^{x) 



dP 

Wix) = -va{x)^, (12) 



where W(x) is the same heating rate as in Eq. IH Then particles are introduced at xa 
with the temperature corresponding to Pq^xa) and allowed to propagate, get injected and 
accelerated. 

This procedure requires a differentiation between thermal and superthermal particles in 
order to calculate the total momentum and energy fluxes upstream (Eqs. Eland [9]). Having 
placed the shock at x = (with the upstream at x > 0), the definition adopted in the model 
is that a particle is thermal if it has never crossed the coordinate a; = going against the 
flow {i.e., it never got injected), and a CR particle otherwise. Therefore, for x > xa Eqs. [8] 
and |9] are rewritten as 

P{X)U^{X) + Pg{x) + / PxVxfiX, p)d^P = POUI + Pgfi, (13) 

1-P{X)U'{X) + ^^^^ + / KvJix, p)d'p = l-poul + ^^ - Qesc, (M) 

2 7-ly 2 7-1 

where /(x, p) is the distribution function containing only cosmic ray contribution in the 
sense of the above definition. At x > xa, where anisotropics of thermal particle distribution 
are important, Eqs. |8] and [9] are used. 

As one can see, between the location xa and the subshock, thermal particles are only 
heated by compression, but if this compression is rapid, like at the subshock, this compressive 
heating is non-adiabatic. In fact, it turns out that the solutions provided by the model contain 

© 2010 RAS, MNRAS 000.[lti24l 



Comparison of Different Methods for NLDSA 11 

a region of finite thickness > x > Xcrit, in which this non-adiabatic compressive heating 
takes place: this region is interpreted as the subshock. 



4 SEMI- ANALYTIC APPROACH 

The semi- analytic formahsm for the NLDSA problem developed by Caprioli, Amato and 
Blasi (hereafter CAB) couples the conservation of mass, momentum and energy flux with 
an analytical solution of the sta tionary diffusion-convection eq uation for the isotropic part 



of the CR distribution function ( ICaprioli. Amato. Blasi 



2010bf l. 



More precisely, the transport equation is taken as in Eq. [5l with the velocity of the 
scattering centres u^ neglected with respect to the fluid one, their ratio being typically of 
order oi va/u <^ 1: 



u 



dl 
dx 



d_ 
dx 



D{x,p) 



dx 



pdudf 

+ o3-7^+^^^'^^ 
3 ax op 



The injection term Q{x,p) is written following iBlasi. Gabici fc Vannonil ( l2005l ) as 

VnoUo 



2 -S{p - Pinj)S{x) 



(15) 



(16) 



Pinj 

where r] is the fraction of particles crossing the shock and injected in the acceleration process 
and Pinj is the injection momentum (see §4.11 for further details). 

Eq. [T5]is solved imposing the upstream boundary condition f{xo,p) = as in KJ model, 
which mimics the presence of a free-escape boundary at a distance Xq upstream of the shock. 
A net flux of particles (pescip) escaping the system from this boundary is expected to grant t he 



achievement of a stationary configuration, as discussed in 



Caprioh. Blasi. Amatol f l2009br ) 



An excellent approximate solution of Eq. [15] is found to be: 



f{x,p) 

4>esc{p) 

where 
A{x,p) 



fship) exp 
'df 







U[X 



D{x',p} 



1 - 



A(x,p) 

Mp) 



-D{p) 



dx 



J Xo 



Up f ship) 

' Mp) 



dx'- 



Uq 



exp 







ji\ 



U[X 



da;" 7^7 r 

D{x",p) 



Ao{p) = A{xo,p) 



(17) 
(18) 

(19) 



D{x',p) 

embeds all the information on the escape flux and the CR distribution function at the shock 
position fship) can be implicitly written as 



fship) 



Wo 



3R 



■tot 



^T^Plij RtotUpip) - 1 



exp 



3R 



tot 



p dp^ 

,,„^. y Rtotu.ip') - 1 



u,{p') - 



4>e 



Uo fship') 



,(20) 
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12 Caprioli, Kang, Vladimirov and Jones 
Rsub 1 /" , „ , X dU 



Up{p) = -jf^-TT-^ dx/(x,p)— . (21) 

Rtot fsh{p) Jxo dx 

Here we introduced the normalized fluid velocity U = u/uq and the compression ratios at 
the subshock Rsub = Ui/u2 and between xq and downstream Rtot = uq/u2. 

The solution of the problem is found using a recursive method, starting from a guess for 
f/i = Rsub/ Rtot (the compression ratios are linked together by Eq. [28]), along with a guess 
fo r the CR distribution function at the shock and the escape flux (the test-particle solution 



of 



Caprioh. Blasi. Amatol l2009bl ). 



• At each step of iteration the pressure in CRs is calculated from the distribution function 
as: 

Pc{x) = — / dpp^'vip) f{x,p), (22) 



in J 



and it is normalized to the valued of Pc^i deduced from the momentum conservation equation 
integrated between positions immediately upstream of the subshock and Xq] namely 

^ = l + ^-C'.-n,.(C^O. (23) 

where Pg^i{Ui) is given by Eq. |271 

• With this normalized CR pressure, it is possible to calculate at first U{x) from momen- 
tum conservation and eventually a new f{x,p) via Eqs. [T7142T1 

The steps above are iterated, leaving f/i unchanged until convergence is reached, i.e. until 
the factor normalizing Pc,i does not change between one step and the following: in general 
this factor is different from 1, which means that the chosen Ui does not provide a physical 
solution of the problem. The whole process is thus restarted with a new choice of f/i until the 
iteration on f{x,p) leads to a distribution function which does not need to be re-normalized 
in order to satisfy momentum conservation: the f{x,p) and the U{x) found in this way 
represent the solution of the NLDSA problem. 

Once the solution has been found, the equation for the conservation of the energy is 
satisfied as well, and the energy carried away by escaping particles, normalized to PqUq/2, 
can be written as 

fi„^rd /''''^""')f"M . (24) 

or also as the difference between the gas energy flux at Xq and the whole energy flux down- 
stream, namely: 
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(25) 



Rl^t M^{-f - 1) Rtot 7-1 Poul Rtot 7c - 1 Poul ' 
where jc — 4/3 is the adiabatic index of the relativistic gas of cosmic rays. 

The time required to run this iterative scheme is less than one minute on a current 
generation CPU. 



4.1 Injection mechanism 

Also in this approach the injection of particles into the acceleration mechanism is taken as 
occurring immediately dow nstream of the shock due to th ermal leakage, as in KJ. How- 
ever, the implementation by iBlasi. Gabici fc Vannonil ( l2005l ) adopted here is simplified with 
respect to the KJ model, because it does not include a smooth transparency function. 

The injection momentum pinj is in fact chosen as a multiple of the thermal momentum of 



downstream particles, namely: Pmj = i^inj—U2/c)pth,2, withpi/j = ^j2mpkBT2. The shift U2/C 
comes from the assumption that thermal particles downstream have a Maxwellian spectrum 
in the fluid reference frame, while Pinj and the equations above are written in the shock 
frame. Under this hypothesis the fraction rj (Eq. [T6l) is uniquely determined by the choice 

oiiinj, namely: 
4 



30F 



(/?,,, -i)ev 



(26) 



4.2 Turbulent heating 

As in the previous approaches, precursor h eating due to the dissip ation of Alfvenic turbulence 



can be taken into account as described in 
tion proposed by 



Amato fc Blasil (120061 ). following the implementa- 



Berezhko fc EUisonl (Il999[ ). In particular, the entropy conservation for the 
gas (Eq. [T2I) determines the gas pressure in the presence of a non-adiabatic heating as 

•l_f/7+l/2(a.)- 



Pg{x) - P.^Ux) [1 + H{x)] ■ H{x) = 7(7 - 1) 



M'o 



M 



A,0 



7 + 1/2 



(27) 



where Pg^ad{x) / {pouD = U^'^{x)/{jMq) is the plasma pressure as calculated taking into 
account only adiabatic compression in the precursor, i.e. Pg^ad{x)/p'^{x) =const. This ex- 
pression serves as an equation of state for the gas in the presence of effective Al fven heatiri g 
and leads to the following relation between the compression ratios Rgub and Rtot (jBlasi 

7 + 1 - Rsubil - 1 



2002): 



p7+l ^^^0-^sub 

R ~ 



Hot 



;i+^i) 



(28) 



and in turn 
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7 + 1 
^*°* = 2f/7(l + i7i)/M2 + ^1(7-1)' ^^"' = ^^^*°* ^^^^ 

where Hi is given by eq. [27] evaluated upstream of the subshock. 

For completeness, the downstream temperature is thus calculated to be 

T. = To (1^) '" -' + ;-'"- ]\f^ (1 + //O . (30) 

\RsubJ 7 + 1 - (7 - l)/?s«b 



5 A BENCHMARK CASE 

In this section the three approaches described above are applied to the same physical prob- 
lem, in order to emphasize analogies and differences among themselves. The benchmark case 
studied here is a plane, parallel shock with velocity Uq = 5000 km s~^ propagating in a ho- 
mogeneous background with magnetic field Bq = 3/iG, particle density po/nip = 0.003 cm~^ 
and temperature Tq = 2.02 x 10^ K, corresponding to a sonic Mach number Mq = 30. The 
corresponding Alfven speed is, va ~ 120 km s~^, with an Alfvenic Mach number. Ma ~ 42. 
The upstream free escape boundary is located at xq = 3.13 x lO^^cm and the diffusion is 
taken as Bohm-like diffusion in the background magnetic field Bq. In approaches where a 
diffusion coefficient is needed, this is equivalent to 

Dip) = D,-^ , (31) 

nipC 

with D^, = 1.043 X 10^^cm^s~^(-Bo/3/iG)^^. In order to give an idea of the length-scale 

imposed here, the choice of Xq corresponds to the diffusion length A = D{p)/uq of a particle 

with momentum p ~ lO'^GeV/c, which approximates the upper cutoff momentum, Pmax, 

obtained in all the solutions. 

This set of environmental parameters is well defined for the stationary approaches of EV 

and CAB, but is slightly trickier for the time-dependent case of KJ. As noted previously, their 

simulations start with a purely gasdynamic shock established by R-H conditions to give a 

shock speed, uq = 5600 km s~^ at rest at x = 0. There are no pre-existing CRs, i.e., Pc{x) = 

at t = 0. As suprathermal particles are injected at the subshock and accelerated, the CR 

pressure increases, while total compression increases and the shock is slowed down. When 

the shock has reached the self-consistent dynamical equilibrium state, the speed becomes 

Mo ~ 5000 km s~^, which equals to the shock speed of steady state solutions obtained by 

the other two methods. We emphasize again that a gasdynamic shock evolved into a CR 

modified shock can reach a true steady state only when CRs escape the shock structure and 
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Figure 1. Particle distribution /sh(p/mc) at the shock position, in the shock frame, multiphed by {p/mc)'^, calculated within 
different approaches as in the legend. 



carry away energy at a rate that balances the rate at which CRs gain energy through DSA. 
In these simulations that escape conies through imposition of a free escape boundary at Xq. 
In EV approach, there is no distinction between thermal and non-thermal particles, 
hence particle injection is intrinsically defined by the choice of the scattering properties (as 
in §3.2p . and so it is not controlled with a free parameter. In order to match this intrinsic 
recipe for the injection embedded in the Monte Carlo approach, we chose values of e^ and 
,^, which lead to similar total compression ratios as EV's, Rtot = P2/ Po- Quantitatively 
speaking, KJ chose e^ = 0.28 (see §2.31) . while CAB fixed ^ = 3.1 ( §4.ip . corresponding to 
the injection of a fraction r^ ~ 3 x 10^^ of the flux of particles crossing the subshock. It turns 
out these choices push the CR acceleration close to its maximum efficiency, because the non- 
linear feedback and shock modification saturate with respect to increased injection. It has 
been shown previously that, for moderately strong shocks (M > 30), the CR acceleration 
efficiency and the shock modification become insensit ive to the values of en and 1^ for the 
injection fraction above a critical value, rj > 10'^ (e.g. jKang fc Joned 120071 ). 
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5.1 Particle spectrum 

In Fig. [1] the particle distribution at the shock position is shown as a function of the normal- 
ized momentum p^, = p/mc, as multiplied by p^ in order to emphasize the main scaling. The 
agreement among the three solutions is indeed very satisfactory: the concave shape of the 
CR spectrum typical of CR modified shock models is recovered within all the approaches: 
in particular, f{p) is steeper than p~^ at low momenta and then flattens towards higher 
momenta, until it turns down sharply above Pmax- 

The peaks of the distributions of the thermal particles calculated with the three differ- 
ent approaches match very well, and the thermal distributions in the two thermal leakage 
injection models are in excellent agreement. In the Monte Carlo approach the transition 
between thermal and non-thermal particles is much broader than in thermal leakage ones. 
Moreover, the transparency function implemented in KJ thermal leakage modeling makes 
this transition smoother than in CAB one, which shows a sharp break between thermal 
and non-thermal populations. Nevertheless, despite of these different ways of dealing with 
particle injection, above 0.1 GeV/c all spectra are very close to each other. 

As estimated at the beginning of the section, the spectrum has to be cut-off around 
Pmax ~ lO^GeV/c as a consequence of the escape of the highest energy particles, which are 
not able to diffuse back to the shock. The position and the shape of the momentum cut-off 
is consistent among the different models, with the EV solution extending about 10 per cent 
higher in energy with respect to the other cases (see also the left panel of Fig. [2]). This is due 
to a subtle difference in the definition of the transport properties of particles: the models of 
CBA and KJ use a diffusion coefficient to describe relativistic particle propagation, and the 
model of EV defines a small-angle scattering rate instead, which leads to some 20 per cent 
difference in the effective particle mean free path, which is a s mall factor compa. red to the 



uncertainty in the model of particle diffusion (see e.g. § 3.1.2 in IVladimirovll2009[ ) 



5.2 Escape flux: spectrum and anisotropy 

The spectrum of particles escaping the system from the upstream boundary (pesdp) is shown 
in the left panel of Fig. [2] (heavy lines), as multiplied by pt/uo and compared with the highest 
end of the spectrum at the shock location, ptfsh{p) (thin lines). The spectra obtained with 
different methods are in good agreement, with a slight shift towards higher energies of the 
peak of the spectral distribution found within the Monte Carlo approach as discussed above. 
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This fact is of great interest since particles escaping the system from the upstream are 
expected to be re sponsible for the highest ener gy region of the galactic diffuse cosmic ray 



spectrum (see e.g. 



Caprioli. Aniato. Blasi.2010a) . In particular, provided that escape occurs 



when the shock evolution is quasi-stationary (a reasonable assumption during the Sedov- 
Taylor stage of a SNRs, when the instantaneous maximum energy is thought to be limited 
by the accelerator size and not by the acceleration time), the rate and the spectrum of 
accelerated particles that are injected in the Galaxy as CRs is found to be the same in all 
the approaches. 

This remarkable agreement between models that solve the particle transport via Monte 
Carlo techniques or via the diffusion equation is not trivial at all, since the anisotropy of 
the particle distribution function is treated in radically different ways. In order to better 
understand this point we introduce the quantity 

J-1 d/xJo 47rpMp/(x,p,/i) 
which measures the anisotropy of the CR distribution function, averaged over p, at the 
position X. Here /i = Px/Pi with /i = +1(— 1) corresponding to particles moving against (in 
the direction of) the flow, and g{x,fi) is normalized such that /_^ g{x,fi)dfi = 1. 

While the Monte Carlo approach directly takes the anisotropy of the distribution function 
into account (see §3.2p . the other approaches solve the diffusion-convection equation for the 
isotropic part only. More precisely, in the diffusive approximation the distribution function 
is taken as the sum of an isotropic term /'■^^(p) and a dipole-like anisotropic one ^f^^^p), 
in such a way that 

f{x,p)c^f\x,p)+fif\x,p) (33) 



This is expected to be a very good approximation for non-re 



ativistic shocks since the n 



th 



SkiUine 



19751 . for a thorough 



term in the anisotropy expansion is of order (u/c)"' (see e.g. 

derivation of the diffusion-convection equation). The introduction of the spatial diffusion 

coefficient in fact allows one to eliminate f^{x,p), and the diffusive flux reads 

- D{p)^^^^^ = f_'^ /(P)c^d/i = \cf\x,p). (34) 

We are particularly interested in the anisotropy at the free-escape boundary, where it is 
expected to be the largest, but there is a caveat in carrying out this analysis. Whichever is 
the value of the distribution function at Xq as a function of p, the anisotropy immediately 
downstream of Xq, which we will refer to as g{xQ,^), and immediately upstream of Xq, 
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Figure 2. Left panel: flux of particles escaping tlie system from the free-escaping boundary at xq {(peso, heavy lines) and 
particle distribution at the shock (/s^, thin lines). Different curves correspond to different models, as in the legend of Fig. [T] 
Right panel: angular distribution, averaged over momenta, of escaping particles (gesc(M)i top) and of the particles immediately 
downstream of the free-escape boundary (g{xo,fj,), bottom). 

indicated with gesc{p,f^), are intrinsically different. This must be true, since upstream of xq 
particles must escape the system, and hence, QesciVilA has to vanish for jj, < uq/c <^ 1. 
Downstream of Xq, instead, particles can diffuse back to the shock and thus Eq. [5^ safely 
holds also for yU < 0. 

Keeping this difference in mind, we can write the anisotropy of the distribution function 
at the position xq also for diffusive approaches using Eqs. [32llMt obtaining 



^(a^o,^) = 2 (1 + ^/^) ' ^ 



/o°° dp p^f\xo, p) _ 3 J^ dp p^(f)esc{p) 



(35) 



(36) 



/o dpp^p{xo,p) 2c /q dpp^p{xo,p)' 
while the anisotropy of the escaping particles reads 

where ^{fi) is the usual Heaviside step function. 

This procedure allows one to work out the expected anisotropy close to the escape bound- 
ary within the numerical approach of KJ, where the boundary condition f{x > xo,p) = is 
imposed and thus f{xo,p) is finite, but also in the semi-analytical approach of CAB, where 
f{xo,p) is imposed to be exactly 0, simply by taking the left limit for x — )■ Xq in the equations 
above. 

The results are shown in the right panel of Fig. [21 along with the result of the Monte 
Carlo simulation by EV at xq (top panel) and at x = 0.99xo (bottom panel). Also in this 
case the agreement among the different approaches is remarkable and, most interestingly, 
the degree of anisotropy which the methods based on the solution of the diffusion-convection 
equation are able to retain is indeed very accurate, even when compared to the Monte Carlo 
method which deals with the fully anisotropic distribution function and which can track the 
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Figure 3. Spatial dependence of the fluid velocity (normalized to mq, 'e/i panel) and density (normalized to pgright panel) in 
the three difl'erent approaches. Here and in the following figures the shock faces to the right and is rest at a; = 0; the upstream 
is truncated at x = 10~*a:o and the downstream is assumed as homogeneous and not in scale. 

smooth transition from the hnear to the step-hke functional form of g{x,n) in the range 
0.99 < x/xo < 1. 

Quantitatively speaking, we can safely conclude that, for plane non-relativistic shocks, 
the discrepancies in the spectra of accelerated and escaping particles between time- dependent 
and stationary approaches or between NLDSA models taking or not into account the anisotropy 
of the distribution function are well within any uncertainty induced, for instance, by any 
observational constraint of the several environmental parameters of a given astrophysical 
object. 



5.3 Hydrodynamics 

The spatial structures of the modified shocks are quite consistent among the three models. 
This is evident from Fig. [3] which compares velocity and density profiles. As stated above, 
the CR acceleration efficiencies are very close to the upper limits allowed by the system when 
turbulent heating is taken into account. The slightly different saturation levels achieved in 
different approaches might thus be regarded as due to intrinsic properties of the models. 

All the models return the same ratio of unperturbed to downstream fluid density Rtot = 
P2/P0 within a tolerance of less than 10 per cent. In particular, the total compression ratios 
are found to be Rtot,cAB — 7.2, Rtot,Kj — 7.3 and Rtot,EW — 7.6, corresponding to CR 
acceleration efficiencies of about 60 per cent of the bulk pressure for CAB and KJ approaches 
and of about 75 per cent in the EV case. The higher shock modification in the EV case is 
due to a larger amount of energy in supra-thermal, non relativistic particles compared to 
the CAB and KJ approaches (see Fig. [1]). 
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Figure 4. Spatial dependence of the gas pressure {left panel) and temperature {right panel) in the three different approaches. 
The spatial coordinate is taken as in Fig. [3] 

In the precursor the total spread among the model profiles is within a 10-20 per cent 
tolerance. Unfortunately, at the moment there are no observations able to resolve the spatial 
structure of a precursor. In addition, all the hydrodynamical evidence of efficient CR accel- 
eration at SNR shocks comes from estimates of th e separation distance between the the for- 



ward shock and the contact discontinuity (see e.g. 



2008 



Warren et al. 



2005 



Cassam-Chenai et al. 



, for Tycho and SN1006, res pectively), or by an estimated supp ression of the down- 



stream plasma temperature (e.g. 



Decourchelle. Ellison fc Ballet 



20001 ) with respect to the 



pure gaseous case. On the other hand, it was suggested th at the Ha profi 



Wagner et al 



es in Tycho's SNR 



20091 ) 



shock might indicate the presence of a CR precursor (e.g. 

In Fig. m the gas pressure and temperature are shown for the three models, confirm- 
ing once again that all the approaches lead to consistent results. The lower values of the 
downstream gas pressure and temperature in the Monte Carlo approach with respect to the 
other cases may be due to the slightly more marked shock modification, since the larger is 
the amount of energy channelled into accelerated particles, the lower is the energy left for 
heating up the downstream plasma. On the other hand, above the thermal peak the Monte 
Carlo solution shows significant deviation from a pure Maxwellian distribution, because the 
injection model adopted in the simulation and the nature of the calculation allow mildly 
accelerated particles (those having crossed the shock only a few times) to be present in the 
shock vicinity. 

Finally, the fraction of the bulk energy density flux pqu^/2 carried away by particles 
escaping the system at the upstream boundary is found to be -Fe^cKj = 0.29, F^scFN = 0.27 
and -FescCAB = 0.23, consistently with the small differences in the acceleration efficiencies 
outlined above. 
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6 COMMENTS AND CONCLUSIONS 

In this work we conducted a comparison among numerical, Monte Carlo and semi-analytical 
approaches to NLDSA theory by solving a benchmark case for an astrophysical plane, non- 
relativistic strong shock. In particular we examined a case in which cosmic rays are efficiently 
accelerated and strongly modify the shock with their dynamical back-reaction. 

The results obtained with the three different techniques that assume specific and dis- 
tinct recipes for physical ingredients such as particle injection, particle escape and diffusion 
properties are very consistent with each other in terms of both accelerated (and escaping) 
spectra and hydrodynamics. These findings represent an actual investigation of the role that 
specific pieces of information embedded in each approach have in modeling cosmic ray mod- 
ified shocks. In particular, we showed how the solution obtained when the stationary state 
is achieved by (a) letting the system evolve at the presence of a free-escape boundary (the 
numerical approach, see^, (b) a steady-state solution in which the particle transport is 
described statistically through a Monte Carlo technique (§3]), or (c) the semi-analytic solu- 
tion of the stationary diffusion-convection equation for the isotropic part of the cosmic ray 
distribution function (§!]) show an agreement well within any uncertainty intrinsic to the 
basic physics of the process. Thus, the agreement would not allow distinction of the methods 
using current observations. 

At the same time, it is worth recalling that all the models are rather flexible and may be 
updated with any new physical insight which may come from present and future observations. 
Very generally, these methods are limited more by our real knowledge of the mechanisms 
going on in a SNR shock than by the possibility of implementation of these effects. 

An important example of a piece of information that can be embedded in all the methods 
is a more detailed description of the magnetic field responsible for the particle scattering. 
In the calculations above we did not include any process leadi ng to the amplif ication of 



the background magnetic field due to CR streaming instabilities (JBell 



dynamical feedback of the magnetic turbulence in the shock dynamics (jCaprioli et al.ll2008l ). 



1978 



20041). nor an v 



On the other hand, by assuming Bohm-like diffusion we have intrinsicall y taken into account 



a mag netic configuration quite different from the typical interstellar one (jLagage &: Cesarsky 



1983a 



b|). Moreover, we included a phenomenological model for a non-adiabatic heating of 



the upstream background plasma as due to the instantaneous damping of all the magnetic 



turbulence produced via resonant streaming instability, as described in 
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(119821 ) ■ This very efficient damping, while largely adopted in order to prevent dramatic 
shock modification, is however inconsistent with magne tic fields of order 1 00 — 500/iG in- 
ferred in the downstream of some young SNRs (see e.g. IParizot et al.l 120061 ): in these envi- 
ronments the shock modification may be controlled by the non-negligible dynamical role of 
the amplified magnetic field, an eff e ct which can be easily included in the models above (see 



Vladimirov. Ellison fc Bykov 



20061 : ICaprioh et al.ll2009ah . 



Another interesting finding of our comparison is that not only Monte Carlo simulations, 
but also methods based on the solution of the diffusion-convection equation can retain 
reliable information about the anisotropy of the distribution function, even close to the free- 
escape boundary, in many astrophysical cases. This fact allows one to use the so-derived 
diffusive CR current to calculate the excitation of modes relevant for the growth of the 
magnetic turbulence, as predicted in many scenarios involving the self-generation of an 



amplified magnetic fie 



Reville. Kirk k Duffvl 



d res ponsible in which particles are scattered very efficiently (see e.g. 



2009 



for the calculation of the Bell non-resonant streaming instability 



excited by escaping particles). 

In this work we considered a situation in which injection occurs very efficiently, as sug- 
gested by the Monte Carlo modeling of the subshock. This cho ice, while on one hand perfectly 
repro duces the observations of the Earth bow shock (e.g. 



Ellison. Mobius &: Paschmann 



1990l ). on the other hand is not firmly supported by observational evidence in SNR shocks. 



For typical values of temperatures and magnetic fields in SNR environments, in fact, the 
region of the spectrum between the thermal and the non-thermal distributions does not pro- 
vide, at the moment, any strong signature testable with multi-wavelength studies of SNRs. 
In this respect, a simple parametrization of the injection process like the thermal leakage 
captures most of the (non-linear) physics which is expected to be relevant in NLDSA theory. 
Further improvements in discrimination between the various models of particle injection can 
be expected from the ongoing work on particle-in-cell simulations of non-relativistic shocks. 

Dealing with strongly modified shocks is indeed very challenging for fully numerical 
methods. When huge gradients are found in the precursor, the system very likely devel- 
ops some instabilities, whose nature may be not only numerical, but also physical. For in- 



stance, the acoustic i nstability studied by lDrury fc Falld ( 1l986l ): lKang. Jones fc Ryul (119971 ): 



Wagner et al 



(|2009[ ) can develop strong inhomogeneities, which can drive turbulence in 
multi-dimensional flows. Numerical methods are thus tools necessary in order to unravel 
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these kind of questions, being in principle able to follow the growth of instabilities even in 
the non-linear regime, typically inaccessible to analytical approaches. 

A numerical simulation is currently the only way to deal with time- dependent NLDSA 
problems. We recall that there are two main ways of accounting for the dependence on time 
of the physical quantities: on one hand we have cases in which the shock velocity is almost 
constant and only CR-related quantities are expected to vary on time scales comparable 
with the acceleration time at Pmax-, and on the other hand we have cases in which also 
"environmental" quantities, like the shock velocity, change with time. 

The former case is realized for instance during the ejecta dominated stage of a SNR, 
in which the shock is not dramatically slowed down by the inertia of the swept-up mate- 
rial. Nevertheless, if a spatial boundary condition is imposed in such a way that particles 
are free to escape from an upstream boundary, also in this case a steady-state solution 
consistent with stationary approaches may be recovered. Otherwise, if Pmax is free to grow 
unconstr ained by spatial limits, the system evolves accord ing to a self-similar behavior as 
shown by iKang fc JonesI (120071 ) ; iKang. Ryu fc Joned ([20091) : also in this case stationary ap 
proaches, with a proper age- limited maximum momentum (JBlasi. Amato fc Capriolil 12007 
see), should be able to reproduce the correct results. 

The latter case, instead, i s expected to be realized in the Sedov- Taylor stage of the evo- 
lution of SNRs. As shown by lCaprioli. Blasi. Amatol (l2009bl ). during this stage the decrease 
of the shock velocity is expected to lead to a decrease of the magnetic field amplification 
and in turn to a drop of the confinement power of the system. At any given time particles 
with momentum close to Pmax{t) cannot diffuse back to the shock and escape the system 
upstream, thus, carrying away a sizable fraction of the energy flow into the shock. This 
situation is thought to lead to the achievement of a quasi- stationary configuration in which 
particle acceleration and losses balance themselves; also in this case a time-independent 
approach might be a good description of the problem. 

In the case in which the assumption of quasi-stationarity is satisfied, modeling an astro- 
physical object might require one to run the NLDSA code many times: one for each time-step, 
for any given choice of the many environmental and model parameters. It is clear that such 
an investigation can be conducted only by using very fast solutions of the NLDSA problem, 



and se mi-analytical approaches can do this job very well (see e.g. 



Caprioli. Amato. Blasi 



2010al . for an example). 



Finally, Monte Carlo codes allow one to take into account strong anisotropics of the 
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CR distribution function, for which the diffusive approximation and hence the solution of 
the diffusion-convection equation for the isotropic part of -F(x, p) only may fail by a large 
amount. This poi nt is of central interest especially for the study of mildly-relativistic and 



relativistic shocks ( jEllison fc Doublell2002l ). since anisotropic corrections to Eq. |5]are of order 
0{u^/c^), as outlined in J JOl 

The last consideration is about an assumption that is ubiquitous in NLDSA theories; 
namely, that particle are efficiently scattered by magnetic irregularities and in a diffusive 
way. This assumption may be broken close to xq, beyond which particles are expected to 
stream away without diffusing, or even in the precursor if the topology of the niagne t ic fiel d 



2004 . 



20051 ). 



is organized in structures like the cavities predicted by Bell's instability (iBelll 
Nonetheless it is worth stressing that in the literature no non-linear particle acceleration 
methods able to relax the hypothesis that particle are scattered in a diffusive way (which 
is embedded also in Monte Carlo approaches) have been developed yet. Such a problem 
could be in principle addressed by Particle in Cell (PIC) simulations, but the computational 
time required to run a case in which the spectrum of accelerated particles spans many order 
of magnitudes is now, and it is going to be for years, prohibitive even for huge computer 
clusters. 
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